Traveling and resting crystals in active systems 
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A microscopic field theory for crystallization in active systems is proposed which unifies the 
phase-field-crystal model of freezing with the Toner- Tu theory for self-propelled particles. A wealth 
of different active crystalline states are predicted and characterized. In particular, for increasing 
strength of self-propulsion, a transition from a resting crystal to a traveling crystalline state is found 
where the particles migrate collectively while keeping their crystalline order. Our predictions, which 
are verifiable in experiments and in particle-resolved computer simulations, provide a starting point 
for the design of new active materials. 
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Self-propelled particles [1 exhibit fascinating collec- 
tive phenomena like swarming, swirling and laning which 
have been intensely explored by theory, simulation and 
experiment, for recent reviews see [2-4]. In marked con- 
trast to passive particles, self-propelled "active" particles 
have an internal motor of propulsion, dissipate energy 
and are therefore intrinsically in nonequilibrium. Exam- 
ples of active particles include living systems, like bac- 
teria and microbes [5], as well as man-made microswim- 
mers, catalytically driven colloids [6j[7] and granular hop- 
pers |5]. 

If, at high densities, the particle interaction domi- 
nates the propulsion, crystallization in an active system 
is conceivable. It is expected that such "active crys- 
tals" have structural and dynamical properties largely 
different from equilibrium crystals due to the intrinsic 
drive. In fact, there is experimental evidence for active 
crystals, both from observations of hexagonal structures 
for catalytically-driven colloids [9] and honeycomb-like 
textures for flagellated marine bacteria [lOj [11] . More- 
over, recent computer simulations have confirmed crys- 
tallization [T2UT4] and proved that melting of active crys- 
tals differs from its equilibrium counterpart. However, 
though field-theoretical modelling of active systems has 
been widely applied to orientational ordering phenomena 
[2j [15] , there is no such theory for translational ordering 
of active crystals nor has a systematic classification of 
active crystals been achieved. 

Here we present a microscopic field-theoretical ap- 
proach to crystallization in active systems and we pro- 
pose a minimal model which has the necessary ingredi- 
ents for both, crystallization and activity. In doing so, 
we combine the phase-field crystal model of freezing [16] 
with the Toner- Tu model for active systems [17] using the 
concept of dynamical density functional theory p~8j [19] . 
On the one hand, the phase-field-crystal (PFC) model 
as originally introduced by Elder and coworkers p~6j [20] 
describes crystallization of passive particles on micro- 
scopic length and diffusive time scales. When brought 
into connection with dynamical density functional theory 



[2TH24] . the PFC model represents in principle a micro- 
scopic theory for crystallization, and it has been success- 
fully applied to a plethora of solidification phenomena 
[I5|I5D| I25H55]. On the other hand, Toner and Tu [17] in- 
vestigated the onset of collective motion in self-propelled 
systems from a general hydrodynamic point of view. Phe- 
nomenological coupling parameters of this model can in 
principle be justified by dynamical density functional the- 
ory [30], too, but it does not describe crystallization. 

In our PFC model for active systems, we find a wealth 
of different crystallization phenomena. First, we identify 
two different types of active crystals which we call "rest- 
ing" and "traveling" depending on their averaged drift 
velocity. A resting crystal possesses vanishing net par- 
ticle flux whereas a traveling crystal is migrating with 
a nonzero velocity while keeping its periodicity. Start- 
ing from a disordered initial state, a traveling crystal 
is typically formed by a coarse-graining process of do- 
mains. The threshold in the driving strength upon which 
traveling crystals are formed depends on the spontaneous 
local orientational order (as prescribed by the coupling 
parameters of the bare Toner- Tu model): if there is no 
such order, the threshold is finite, while there is no such 
threshold in the presence of spontaneous orientational or- 
der. We further identify a transition from a hexagonal to 
a rhombic traveling crystal if the drive is increased fur- 
ther and find also resting and traveling lamellar phases 
with one-dimensional periodic ordering. Finally the oc- 
currence of honeycomb- like structures can be explained 
as well within our model. The knowledge and control 
over these crystalline states provides an attractive start- 
ing point for the design of novel active materials since 
active crystals possess unique structural, phononic, and 
rheological properties. 

In the following, we first describe our model and then 
discuss numerical and analytical results. Our dynami- 
cal equations are for the local one-particle density field 
ipi (r, t) which is a conserved scalar order parameter 
and basically describes the reduced density modulation 
around a fixed averaged density t/5 as in the traditional 
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PFC model p~6j [20] , and for a polarization vector field 
P(r, t) which describes the local polar ordering. Activity 
enters into the equations via a nonzero self-propulsion 
velocity vq. In suitably scaled units of time, length and 
energy, our basic dynamic equations read 



(i) 



Mi = v 2 |^- Uo v-p, 

Sip! 

d t P = V 2 ^-D r S ^-v ^ 1 . (2) 

Here, D r is the rotational diffusion coefficient of the par- 
ticles and J 7 is a free energy functional of ipi and P gained 
from density functional theory. The equations and 
([2| are consistent with phenomenological symmetry ar- 
guments and involve the simplest nontrivial coupling be- 
tween the two order parameter fields ipi and P. As out- 
lined in the supplemental material [31 , they can also be 
derived from microscopic dynamical density functional 
theory within an appropriate gradient and Taylor expan- 
sion of the order parameter fields [I9j [32] , see also [33] . 
In the sequel, we shall consider two spatial dimensions 
only. 

We now further specify the free energy functional T to 
T = T v f c + F-p where 

T pSc = j d 2 r{\i,[e + {l + V 2 f}^ + \^} (3) 

is the traditional PFC functional [To) l20] describing the 
tendency of the material to form periodic structures. 
Here, e sets the temperature |16[ l20] , and the order pa- 
rameter tp corresponds to the total density tp = tp + xp\ . 
The polarization-dependent part 



ic 4 (P 2 ) 2 } 



(4) 



describes local orientational ordering due to the active 
driving following the approach by Toner and Tu [17] for 
neglected convection. The functional possesses two cou- 
pling parameters C\ and C4 which govern the local ori- 
entational ordering due to the drive. If C\ = C4 = 0, 
only gradients in the density ipi can induce local polar 
order P of the active driving. For C\ > (C4 = 0) dif- 
fusion tends to reduce the polar order generated by the 
density gradients. In the third case, C\ < and C4 > 0, 
a net local driving spontaneously emerges already in the 
absence of density gradients. 

Clearly, on the one hand, for vanishing self-propulsion 
vq = 0, Eqs. ([!]) and Q decouple and the density 
equation reduces to the usual phase field crystal model 
[TBI [20] . On the other hand, if F p f c is neglected, the 
remaining terms are contained in the model by Toner 
et al. [TTJ [34], except for the higher-order term in P 
that contributes to translational diffusion. Summariz- 
ing, Eqs. ([TJj-Q form a minimal approach to characterize 
crystallization in actively driven systems. 
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FIG. 1. Phase diagrams ("rcryst": resting crystals; "rlam": 
resting lamellae; "tcryst": traveling crystals; "tlam": trav- 
eling lamellae), (a) For C\ — 0, vo = the equilibrium 
phase field crystal model is recovered. Equilibrium phase 
boundaries given by the energy functional are indicated for 
the liquid-hexagonal (dashed line) and hexagonal-lamellar 
(dash-dotted line) transitions, (b) For C\ — 0.2, vo = 0.35 
the structures are still at rest, but the phase boundaries are 
shifted by a value As. (c) For C\ — 0.2, vo = 0.7 the struc- 
tures are traveling and phase boundary lines are omitted for 
clarity. The black stars in the bottom left of panels (b) and 
(c) mark the intersection points with the curve in Fig. [3] In 
all cases C4 = 0, D r = 0.5. 



We numerically determined the phase diagram by scan- 
ning the ijj-e plane while keeping the parameters Ci, C4, 
and fixed. As for any numerical result reported sub- 
sequently, we proceeded in the following way. For each 
set of parameter values £, Ci, C4, i>o) we started from 
random initial conditions and then iterated Eqs. Q-Q 
forward in time. Numerical measurements were carried 
out after equilibration, and a systematic finite size study 
was performed to test the validity of our results. 

For the decoupled case = 0, the equilibrium phase 
diagram [20 corresponding to the energy functional 
Eq. ([3| is shown in Fig. [TJa). For nonzero active drive 
^o, we will first report on the case C\ > 0. 

When we moderately increase from zero for C\ > 0, 
the phase boundaries undergo a temperature shift As to 
lower temperatures. An example is depicted in Fig.JTJb). 
Comparison to Fig. [TJa) shows that switching on the ac- 
tive drive melts crystals and lamellae close to the liquid 
phase boundary. The patterns still remain at rest, how- 
ever. For this case, a linear stability analysis and derived 
amplitude equations for ipi and P predict Ae oc Vq/C±, 
which was also verified numerically. 

We present an example snapshot of the resting crys- 
talline phase in Fig. [2|a). The peaks of the density dis- 
tribution ipi form a hexagonal lattice as dictated by the 
PFC energy functional. P points down the density gra- 
dients. Consequently the polarization field forms "+1"- 
defects centered at the density peaks. Since P describes 
the local direction of active drive, density is convected 
out of the peaks by the active propulsion vq. This mech- 
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FIG. 2. Snapshots of the order parameter fields for different 
phases that are observed when increasing the active drive vo 
at (^,£,Ci,C 4 ) = (-0.4,-0.98,0.2,0): (a) resting hexagonal, 
vo = 0.1, (b) traveling hexagonal, vo = 0.5, (c) traveling 
quadratic, vo = 1, (d) traveling lamellar, vo = 1.9. The 
phases are depicted by plotting the density field ipi : brighter 
color corresponds to higher densities. Thin bright needles 
illustrate the polarization field P that points from the thick 
to the thin ends. In panels (b)-(d) the predominant direction 
of motion is indicated by the bright arrows. Only a fraction 
of the numerical calculation box is shown. 



anism counteracts the density diffusion into the peaks 
described by the PFC energy functional. Therefore lower 
temperatures are necessary for the patterns to form in the 
presence of an active drive, corresponding to the temper- 
ature shift Ae in Fig.JTJb). In the resting crystalline and 
lamellar case, both tendencies balance each other so that 
the averaged net particle flux vanishes. 

When we increase the active drive, we find that the 
density peaks start to travel above a critical value vq^ c . 
Such a state is illustrated in Fig. [5|b). As we can see, 
the centers of the density peaks are now shifted with re- 
spect to centers of the "+1" -defects in the polarization 
field. This reduced symmetry induces active propulsion: 
a net orientation of the polarization field emerges when 
averaged over the area of a single density peak. The con- 
sequence is an active convection of each density peak, 
originating from the PFC density modulation. These re- 
sults are in agreement with a linear stability analysis of 
Eqs. ([!]) and Q which predicts that propagating modes 
appear above a threshold value vo, c - 

With further increasing vq the hexagonal pattern is de- 
formed to a rhombic one. In the end, we observe a nearly 
quadratic structure as depicted in Fig. [5|c). This struc- 
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FIG. 3. Sample- averaged magnitude v m of the crystal peak 
velocities (left scale) and polar order parameter p v of the crys- 
tal peak velocity vectors (right scale) as a function of vo for 
(^,£,Ci,C 4 ) = (-0.4,-0.98,0.2,0). The threshold corre- 
sponds to the onset of collective crystalline motion. Thick 
arrows mark the positions where the snapshots of Fig. [2] were 
taken; the black stars indicate the intersection points with 
the phase diagrams in Figs, [ljb) and[TJc). The region above 
threshold where regular swinging motion could be observed 
is marked in gray. Inset: peak trajectories illustrating a state 
of regular swinging motion in a hexagonal crystal; different 
colors correspond to different peaks; only trajectories of a 
horizontal row of density peaks are shown that started at the 
bottom and were traveling to the top of the picture while 
tracking was performed. 



tural hexagonal-rhombic-quadratic transition appears to 
be smooth and continuous. 

Finally, we observe that the traveling crystal can be 
melted into a traveling lamellar state if vo is increased to 
still higher values. An example snapshot of such travel- 
ing lamellae is shown in Fig. [2^d). In contrast to the 
hexagonal-rhombic-quadratic distortion of the travel- 
ing crystalline lattices, the traveling crystalline-lamellar 
transition occurs rather abruptly. The transition is also 
evident when we compare the two phase diagrams in 
Figs.Qb) andQc). There, with increasing v^, the trav- 
eling lamellar regions grow into the traveling crystalline 
regions. 

To quantify the scenario further, we tracked the motion 
of each density peak. We determined the individual peak 
velocities v$, where i = 1,..,N P and N p denotes the num- 
ber of peaks. Samples of up to 1000 density peaks were 
investigated. The sample- averaged peak velocity mag- 
nitude follows v m = J2^ 1 \\vi\\/N p . In addition, we 
calculated the degree of polar orientational order of the 
normalized peak velocities p v = || (5^=Ti v i/ll v i||) 
This order parameter detects whether the peaks move 
coherently (collectively) into the same direction. 

For C\ > 0, Fig. [3] clearly illustrates the existence of a 
threshold value vq jC at which propagation starts. As indi- 
cated in the inset, we observed a regular swinging motion 
of the hexagonal crystals close to the threshold. With in- 
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FIG. 4. Coarse-graining in a sample of about 1000 travel- 
ing density peaks. The sample is in the traveling hexagonal 
crystalline state. Depicted are instant pieces of all peak tra- 
jectories, drawn in different colors. These pieces connect to 
lines of equal orientation within collectively moving domains. 
Over time, some domains of equal velocity orientation grow on 
the cost of others. Parameter values correspond to Fig. [2jb) . 
Dimensionless times [31]: (a) 4500, (b) 50000, (c) 70000. 



creasing values of active drive, we subsequently find the 
states illustrated in Fig. [2ja)-(d). The averaged peak 
velocity magnitude v m monotonously increases, until it 
abruptly drops at the transition from traveling quadratic 
crystals to traveling lamellae. We can obtain a trav- 
eling quadratic crystal from superimposing perpendicu- 
larly oriented and traveling lamellae. Their intersections 
form peaks that travel a/2 times faster than each single 
lamella by itself, which approximately corresponds to the 
magnitude of the drop in the v m -curve. 

Furthermore, we observe in Fig. [3] that the polar peak 
velocity order parameter p v jumps to a value close to 
one at the threshold and then further increases. This 
indicates that after equilibration of the sample the den- 
sity peaks migrate coherently into the same direction and 
the crystal travels as a single object. However, at each 
value of vq, this collective motion has to first develop 
from the disordered initial state. The latter process oc- 
curs through a coarse-graining dynamics as qualitatively 
illustrated in Fig. [4] The panels depict traveling hexago- 
nal peak trajectories at different times of coarse-graining. 
First, collectively moving crystalline domains form from 
the disordered initial state. However, the migration di- 
rections of different domains are not identical. Over time, 
some domains grow on the cost of others, until a collec- 
tively traveling crystal emerges. 

Finally, if we set C\ < and C4 > 0, a net polar direc- 
tion P of self-propulsion spontaneously occurs as in the 
Toner-Tu model [T7J [34] . For this scenario we never ob- 
served a finite threshold value of vq. Propagating struc- 
tures evolved for all tested nonzero values of vq. Again, 
a transition from hexagonal to rhombic to quadratic and 
then to lamellar structures was observed with increas- 
ing vq. Furthermore, we note that our equations of mo- 
tion are invariant under the transformation ip — — -0, 
ipi — >• —0i, P —> —P. Because of these symmetry re- 
lations, our analysis equally applies for the investigation 
of active honeycomb textures that follow from > 0. 
Such textures were observed for flagellated marine bac- 



teria main]. 

In summary, we extended the phase field crystal model 
[T6| [20] to active systems by combining it with the ap- 
proach of Toner et al. [17] [34]. As a result, the active 
drive favors the liquid and lamellar states in the PFC 
phase diagram and induces a wealth of new active crys- 
talline states of hexagonal, honeycomb, rhombic, and 
quadratic texture. The global motion of all these struc- 
tured states is either "resting" or "traveling" . The transi- 
tion from "resting" to "traveling" involves a complex in- 
termediate swinging motion. When prepared from an ini- 
tially disordered state, traveling crystals emerge through 
coarse-graining from a multidomain texture. 

Our model can be extended from two to three spatial 
dimensions where more crystalline lattice structures be- 
come stable [24] and to binary mixtures of driven and 
undriven particles promising a rich variety of mixed ac- 
tive crystals. In principle, our predictions are verifiable 
in experiments on self-propelled particles and in particle- 
resolved computer simulations at high density [T2HT4] . 
Very recently, traveling crystals have in fact been found 
in such simulations [35]. Since the new traveling crys- 
talline structures show a nontrivial dynamical response, 
they may serve as a building block for a new class of 
active matter with unusual rheological, phononic, and 
possibly also photonic properties. 
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